K-Space Approach in Optical Coherence Tomography: Rigorous Digital Transformation of Arbitrary-Shape Beams, Aberration Elimination and Super-Refocusing beyond Conventional Phase Correction Procedures

For the most popular method of scan formation in Optical Coherence Tomography (OCT) based on plane-parallel scanning of the illuminating beam, we present a compact but rigorous K-space description in which the spectral representation is used to describe both the axial and lateral structure of the illuminating/received OCT signals. Along with the majority of descriptions of OCT-image formation, the discussed approach relies on the basic principle of OCT operation, in which ballistic backscattering of the illuminating light is assumed. This single-scattering assumption is the main limitation, whereas in other aspects, the presented approach is rather general. In particular, it is applicable to arbitrary beam shapes without the need for paraxial approximation or the assumption of Gaussian beams. The main result of this study is the use of the proposed K-space description to analytically derive a filtering function that allows one to digitally transform the initial 3D set of complex-valued OCT data into a desired (target) dataset of a rather general form. An essential feature of the proposed filtering procedures is the utilization of both phase and amplitude transformations, unlike conventionally discussed phase-only transformations. To illustrate the efficiency and generality of the proposed filtering function, the latter is applied to the mutual transformation of non-Gaussian beams and to the digital elimination of arbitrary aberrations at the illuminating/receiving aperture. As another example, in addition to the conventionally discussed digital refocusing enabling depth-independent lateral resolution the same as in the physical focus, we use the derived filtering function to perform digital “super-refocusing.” The latter does not yet overcome the diffraction limit but readily enables lateral resolution several times better than in the initial physical focus.


Introduction
Various methods are used to describe and numerically simulate image formation in optical coherence tomography (OCT), among which much attention has been paid to such a powerful approach as Monte-Carlo (MC), the basics of which were recently discussed in overview [1] and its direct application to OCT can be found, e.g., in [2,3].The strength of MC-methods is related to the very MC principle, in which the image is formed due to a huge number of scattering events of individual photons.This principle implies a rather high computational demand, but at the same time makes the MC approach rather general.In particular, the MC approach allows for the accounting of multiple scattering, which is difficult to introduce in other approaches.
It can be recalled that the principle of OCT-image formation [4][5][6] is based on ballistic single backscattering of the illuminating optical beam (i.e., backscattering described in the first-order Born approximation), so that, among various reasons for OCT-image Sensors 2024, 24, 2931 2 of 20 degradation, multiple scattering is only one such factor, consideration of which for many applications is not of primary importance.In view of this, significant attention was paid to the development of models of OCT-scan formation based on the wave representation and single backscattering of the illuminating beams.In such wave-based approaches, the illuminating field of OCT systems can be described using somewhat different representations.In some cases, important features of OCT-image formation could be understood using simple plane-wave approximations in various forms [7,8].For example, image formation in spectral-domain OCT can be described by representing the illuminating field as a set of co-propagating plane waves with different axial wavenumbers.This approximation, even without special attention to the lateral beam structure, can be sufficient to describe the axial form of the point-spread function (PSF).In some variants, along with the description of the PSF form in the axial direction by considering a set of spectral components with different axial wavenumbers, the lateral form of the OCT beam was also taken into account using the spatial-domain description.Often, the lateral field distribution was described using Gaussian beams [9][10][11][12].Such approaches opened the possibility of fairly accurately describing the lateral form of PSF, including the influence of the illuminating beam focusing.However, such Gaussian-beam-based models do not allow one to describe other beam forms (e.g., Bessel beams), which also represent interest for utilization in OCT [13][14][15].
Alternatively, for describing the lateral distribution of the illuminating and scattered fields, other models based on angular spectra (i.e., representation of the field as a set of plane waves with various directions) have also been developed with one or another degree of rigor (see, e.g., [16][17][18][19][20]).Such approaches combine the spectral representation of both the axial and lateral structures of illuminating beams and, therefore, can be termed a full K-space description.Recently, a comprehensive overview of various K-space-based approaches developed in various techniques of optical visualization (not only OCT in the conventional sense) was published [21].It was pointed out in [21] that the K-space description is a rather general approach, allowing one to analyze very diverse factors in a unified manner, including various beam shapes, effects of dispersion, and aberrations of the optical-wave front.
In application to conventional point-scanning OCT systems, adequate post-processing of the acquired 3D sets of complex-valued OCT signals opens flexible possibilities for various transformations of the initial OCT data.In particular, due to such processing, it becomes possible to introduce corrections for PSF distortions caused by aberrations and/or dispersion effects during wave propagation [21].It also becomes possible to perform digital refocusing by reassembling the spectral components of the illuminating field with corrected phases, such that the depth-independent focusing can be formed over the entire 3D dataset (see examples of such digital transformations in [19,21]).Similarly, the problem of digital corrections of aberrations was discussed (e.g., [21,22]), so that the quality of images obtained by an OCT system with an inexpensive optical system and appreciable aberrations can be made comparable with the quality of images obtained with a much more expensive aberration-less optics.
Here, using the terminology of paper [21], we consider a modified version of the K-space description of OCT-image formation in the development of recent paper [19].The main limitation of this approach is the use of single scattering, which is a generally accepted approximation [6,20,21] corresponding to the very principle of OCT operation.We analyze the most common point-scanning scheme of OCT.In this scheme, the illuminating and then backscattered optical signal passes forth-and-back through the tissue and the same illuminating/receiving aperture.In particular, in the presence of aberrations at this aperture, the optical field passes twice through the source of aberrations, so that the aberration-related distortions in the illuminating and scattered fields experience nonlinear mixing at the reception.The forth-and-back propagation of the optical wave is considered in the K-space (angular-spectrum) representation rigorously without such often-used simplifying limitations as the paraxial approximation or utilization of Gaussian beams.Even if the derived integral expressions without such approximations cannot be evaluated Sensors 2024, 24, 2931 3 of 20 analytically, they can readily be estimated numerically for arbitrary shapes of illuminating beams, not necessarily possessing circular symmetry, e.g., in the presence of anisotropic aberrations, as demonstrated below.
The main goal of the present study is the use of the K-space approach, basically similar to [19], to derive a general form of a filtering function that is intended to transform the initial 3D set of complex-valued OCT data into a "desired"(or "target") form.The latter can be understood in a rather general sense.An essential feature of the proposed filtering procedures is the utilization of both phase and amplitude transformations, unlike often discussed procedures based on only phase transformations.After presenting the basic equations describing the initial OCT-signal formation and deriving the general form of the filtering function, we illustrate the generality and efficiency of the proposed approach by considering several instructive examples.In particular, we demonstrate a non-Gaussian example of the transformation of a 3D set of OCT data (initially corresponding to a beam with a Lorentzian profile to the 3D image for a Bessel beam).Examples will be given to demonstrate that the influence of arbitrary (including anisotropic) aberrations at the illuminating/receiving aperture can be digitally compensated (despite the nonlinear mixing of aberrations acquired during forward and backward propagation).Next, in addition to conventionally discussed digital refocusing intended to obtain the same lateral resolution as at the depth of physical focus, we consider an example of digital filtering enabling "super refocusing".This super-refocusing does not yet overcome the diffraction limit imposed by the illuminating-light wavelength, but allows one to overcome the lateral resolution limited by the initial radius of the focal waist of the illuminating beam, although the increased resolution in such super-refocused images is obtained at the expense of some reduction in the signal-to-noise ratio (SNR).

Rigorous K-Space Formulation of OCT-Image Formation in the 1st-Order Born Approximation of the Scattered Signal
The spectral approach in various forms has been widely used for describing OCTimage formation (e.g., [16][17][18]21,23]).In what follows, we will use the K-space formulation presented in [19] as a basis for the further analysis.In the framework of ballistic first-order scattering (which corresponds to the 1st-order Born approximation), this description is rather rigorous and does not need other popular approximations, such as the widely used paraxial approximation or mandatory utilization of Gaussian beams.
Each spectral component of the illuminating beam with wavenumber k is characterized by a complex-amplitude distribution U L (x − x 0 , y − y 0 ; k) at the output lens aperture with the axis passing through lateral coordinates (x 0 , y 0 ).The source is characterized by the discrete spectrum S(k = k n ) over the wavenumber.The lens is orthogonal to z-axis, corresponding to the axial direction, the lens plane has the axial coordinate z = z d as shown in Figure 1.Coordinate z d is the distance between the aperture and the tissue boundary (and often this gap is filled with an immersion layer).Scanning of the beam position with lateral coordinates (x 0 , y 0 ) is used for obtaining the 3D volume of OCT data.We also use a widely utilized approximation of discrete scatterers [24] that are characterized by coordinates (x s , y s , z s ).Such localized scatterers are convenient for representing the initial and transformed images, but the discussed spectral transformations of optical fields are rigorous and independent of the scatterer density.Like in real-spectral-domain OCT, we represent the illuminating beam as a set of spectral components with wavenumbers k = k n .
For a given wavenumber k, the distribution U L (x − x 0 , y − y 0 ; k) of the complex-valued field over the input/output lens with the axis passing through lateral coordinates (x 0 , y 0 ) is characterized by the angular spectrum ÛL (x 0 , y 0 ; k x , k y , k) defined by the following Fouriertransform over lateral coordinates (x, y): where k x , k y are the lateral components of the total wavenumber k; g(k x , k y , k) is the angular spectrum of the illuminating-beam distribution U L (x − x 0 , y − y 0 ; k) over lateral coordinates (x, y) for the axis position (x 0 = 0, y 0 = 0).The representation of the illuminating beam via its angular spectrum g(k x , k y , k) has the advantage that it allows one to describe an arbitrary beam profile as a set of plane waves with the wave vector components  has the advantage that it allows one to describe an arbitrary beam profile as a set of plane waves with the wave vector components ) , , ( During propagation between the input-lens plane and plane s z z = corresponding to the depth of s-th scatterer each of the plane-wave components acquires the additional phase , which means that the complex-valued amplitude of each propagating components is multiplied by the following propagator function Here, subscript "s" means that s h depends on the depth s z of s-th scatterer and it is assumed that ) exp(ikz corresponds to propagation along z-axis (i.e., in-depth propagation).Therefore, according to Equations ( 1) and ( 2 .In contrast, lateral coordinates During propagation between the input-lens plane z = z d and plane z = z s corresponding to the depth of s-th scatterer each of the plane-wave components acquires the additional x − k 2 y , which means that the complex-valued amplitude of each propagating components is multiplied by the following propagator function h s (k x , k y , k): Here, subscript "s" means that h s depends on the depth z s of s-th scatterer and it is assumed that exp(ikz) corresponds to propagation along z-axis (i.e., in-depth propagation).Therefore, according to Equations (1) and (2), the angular spectrum Ûs (x 0 , y 0 , k x , k y , k) of the field propagated from the illuminating aperture with the lateral coordinates (x 0 , y 0 ) to the location of the s-th scatterer with coordinates (x s , y s , z s ) is given the following product Notice that spectrum Ûs (x 0 , y 0 , k x , k y , k) in Equation (3) depends only on the axial coordinate of the scatterer via function h s (k x , k y , k).In contrast, lateral coordinates (x s , y s ) of s-th scatterer do not enter the incident-wave spectrum (3).Its inverse Fourier transform over spectral coordinates (k x , k y ) gives the complex amplitude U s (x 0 , y 0 ; k) of the incident beam with lateral coordinates (x 0 , y 0 ) in the location (x s , y s , z s ) of the scatterer: x +k 2 y ≤k 2 Ûs (x 0 , y 0 ; k x , k y , k)e i(k x x s +k y y s ) dk x dk y = = We emphasize that subscript "s" in U s (x 0 , y 0 ; k) means that the incident-wave amplitude depends not only on the position of the beam axis (x 0 , y 0 ), but also on coordinates (x s , y s , z s ) of the s-th scatterer.
Propagation of the backscattered signal can be conveniently represented, as shown in Figure 2, in which the back-propagation path is mirrored, as if the back-propagated signal effectively propagates in the same direction along the Z-axis.In this represen-tation, the symmetry between forward and back propagation is especially clear.The distance of the two paths is exactly the same, and the axial back-propagation from the s-th scatterer towards the lens again corresponds to multiplication by the propagator h s (k x , k y , k) given by Equation (2).Next, for each position (x 0 , y 0 ) of the lens axis, when passing through the lens, the amplitudes of received spectral components are multiplied by the factor g(k x , k y , k)e −ix 0 k x −iy 0 k y given by Equation (1).Thus, the spectrum transformation during the back propagation is described by the same total spectral factor h s (k x , k y , k)g(k x , k y , k)e −ix 0 k x −iy 0 k y as in Equation (3) for the forward propagation.Consequently, the transition from the spectral representation to the received signal amplitude is again given by the Fourier transform, similar to Equation (4).
plitude depends not only on the position of the beam axis ) , ( 0 0 y x , but also on coordinates ) , , ( of the s-th scatterer.Propagation of the backscattered signal can be conveniently represented, as shown in Figure 2, in which the back-propagation path is mirrored, as if the back-propagated signal effectively propagates in the same direction along the Z-axis.In this representation, the symmetry between forward and back propagation is especially clear.The distance of the two paths is exactly the same, and the axial back-propagation from the s-th scatterer towards the lens again corresponds to multiplication by the propagator   Notice that such an equivalent scheme should be represented for every particular scatterer depth z s . It should also be taken into account that in the plane z = z s passing through s-th scatterer the spectral components of the scattered signal with the wavenumber k can be represented as K s • U s (x 0 , y 0 ; k), where the scattering coefficient K s characterizes the proportionality between the scattered-wave amplitude and the amplitude U s (x 0 , y 0 ; k) of the wave incident onto the scatterer located at (x s , y s , z s ).Notice that coefficient K s may be different for different scatterers.Taking this remark into account the amplitude of the received signal B s (x 0 , y 0 ; k) from s-th scatterer can be written similarly to Equation (4): so that recalling Equation (4) for U s (x 0 , y 0 ; k), one can rewrite Equation ( 5) in a simple and compact form Here, we assume that for widely used OCT systems with a fairly narrow bandwidth of the illuminating source of the order of several percent and for nearly point-like scatterers with broad scattering diagrams, the dependence of the coefficient K s on the wavenumber k and wavenumber projections (k x , k y ) can the neglected.Consequently, K s can be considered approximately constant.However, if necessary, the above-mentioned effects can be incorporated into the dependence K s = K s (k x , k y , k) and this coefficient should enter the Fourier integral in Equation (5).It should also be emphasized that the symmetry of the received and incident at s-th scatterer signals in Equations ( 5) and (6) holds for beams of arbitrary form, not only Gaussian ones.
To find the total form of an A-scan corresponding to the position of the illuminatingbeam axis (x 0 , y 0 ) the contribution of all scatterers should be summed by performing summation over index s of scatterers, and the (inverse) Fourier transform should be performed to combine the contributions of all received spectral components with the wavenumbers k = k n : Here, the abbreviation DFT (Digital Fourier Transform) emphasizes that this Fourier transform is made numerically, for example, using efficient FFT algorithms realized in many numerical packages; the additionally introduced function S(k n ) describes the amplitudes of spectral components with wavenumbers k = k n of the illuminating-source spectrum.
Equations ( 1)-( 7) constitute a convenient, computationally highly efficient framework for simulating various types of OCT scans for exactly controllable and flexibly variable conditions.Depending on the particular problem, the simulated visualized regions with realistic millimeter-scale sizes may contain a few scatterers for demonstration.Also, configurations with > 10 5 scatterers separated by several micrometers to imitate the spatial density of scatters typical of real biological tissues can be readily simulated, as was recently demonstrated in [19].In particular, simulations of regions with either regularly or stochastically moving scatters can efficiently be made, for which alternative methods (such as Monte-Carlo-based ones) are not well suitable.
Concerning the further discussed procedures for digital transformation of 3D sets of OCT data composed of A-scans given by Equation ( 7), it should be noted that for a given position (x 0 , y 0 ) of the illuminating-beam axis, we obtain only one 1D A-scan A x 0 , y 0 , z q given by Equation (7).The latter contains complex amplitudes B s (x 0 , y 0 ; k n ) of the received wave components with wavenumbers k n .Amplitudes A x 0 , y 0 , z q and B s (x 0 , y 0 ; k n ) received for a given position (x 0 , y 0 ) of the beam axis do not yet give information about the distribution of amplitude and phase of the OCT signal in the lateral directions, which is indispensable if we want to digitally transform the OCT image (e.g., for performing digital refocusing).In the K-space description, such transformations correspond to the reassembling of angular spectral components of the received signals.Therefore, to find the set of spectral components Bs (k x , k y , k n ) = DFT x 0 ,y 0 →k x ,k y [B s (x 0 , y 0 ; k n )] one needs to perform scanning of the illuminating beam in xand y-directions to obtain a sufficiently large and densely spaced set of acquired complex-valued A-scans B s (x 0 , y 0 ; k n ) for various (x 0 , y 0 ).These coordinates should be equidistantly spaced to efficiently perform the Fourier transform from the spatial domain to the K-space.Next, to enable the desired lateral resolution, the initial set of spatial data should be sufficiently dense to satisfy the Nyquist-Kotelnikov criterion.In particular, to enable OCT-image refocusing for a focused beam, the distance between adjacent positions of A-scans in both lateral directions should be smaller than the focal radius of the beam.The corresponding examples will be given in the following sections.

A Remark about Independence of Manipulations with the Angular Spectrum of OCT Signals on the Lateral Positions of Scatterers
In what follows, we perform various transformations (such as refocusing, etc.) with the acquired 3D sets of OCT data.For practical usage, a key point is that, for such transformation of OCT images as refocusing, the procedures certainly should account for the specific depth but should not require a priori knowledge of the lateral positions of scatterers.Mathematically, this means that when performing the corresponding manipulations in K-space with OCT data, the lateral positions of scatterers can be factorized independently of the axial positions of these scatterers.
To demonstrate this, we consider Equation ( 7) describing the complex-valued A-scans in the acquired 3D-pack of OCT data.To manipulate the complex-valued quantities A x 0 , y 0 , z q in K-space, they should be transformed to the spectral form by applying the Fourier transform DFT x 0 ,y 0 →k x ,k y to Equation (7).Since the Fourier transform DFT x 0 ,y 0 →k x ,k y , as well as the summations and the inverse Fourier transform DFT k n →z q in Equation ( 7) are linear operations, their order can be arbitrary.Therefore, in Equation ( 7), we can apply DFT x 0 ,y 0 →k x ,k y to amplitudes B s (x 0 , y 0 ; k) of the received signal from s-th scatterer: Bs (k x , k y , k) = DFT x 0 ,y 0 →k x ,k y {B s (x 0 , y 0 ; k)}.The structure of B s (x 0 , y 0 , k) is given by Equation (6).Bearing in mind that the Fourier transform of a product corresponds to the convolution of the Fourier-transform of the multipliers, one can write: Recalling that U s (x 0 , y 0 ; k) is given by Equation ( 4), one can notice that in Equation ( 8), the convolution can be explicitly represented in the following form where we denote It is clear from the structure of Equations ( 9) and (10) that K s and the lateral coordinates (x s , y s ) of the scatterer enter the Fourier transform of the received amplitude B s (x 0 , y 0 ; k) only as a factor K s e −i[k x x s +k y y s ] independent of the convolution integral I(z s ; k, k x , k y ).The latter is given by Equation (10), in which the axial position z s of s-th scatterer enter via the dependence of function h s (. ..) on z s .

Spectral Form of the Filtering Function to Transform the Initial Set of A-Scans
A(x 0 ,y 0 ,z q ) to Arbitrary Target form A (T) (x 0 , y 0 , z q ) Now, we consider how the above-derived relationships can be used to solve the following problem.Let us assume that the illuminating field has the complex-valued amplitude U L (x, y; k) after the illuminating-lens aperture and the acquired A-scans are A x 0 , y 0 , z q .Our goal is to transform the initial 3D set of A-scans A x 0 , y 0 , z q to another ("target") form A (T) x 0 , y 0 , z q .This means that the target form of scans corresponds to another field distribution U (T) L (x, y; k) at the illuminating aperture: In particular, the distribution U L (x, y; k) may contain some aberrations at the illuminating/receiving aperture, whereas U (T) L (x, y; k) is a desired aberration-free distribution.In other cases, the desired distribution U (T) L (x, y; k) may correspond, for example, to refocusing to another depth.
We recall that the field distribution U L (x, y; k n ) in the K-space is described by the angular spectrum g(k x , k y , k) given by integral (1).Similarly, the target distribution (1).Then, in integral relationship (10) one should substitute g(k x , k y , k) by g (T) (k x , k y , k) to obtain the target convolution function Next, we recall that Equations ( 9) and (10) define the angular spectrum Bs (k x , k y ; k) of the received signal from s-th scatterer for the initial field distribution U L (x, y; k n ).Equations with the same structure as ( 9) and (10) also define angular spectrum B(T) s (k x , k y , k) for the target distribution U (T) L (x, y; k n ).Consequently, the angular spectrum Bs (k x , k y , k) can be transformed to the desired form B(T) s (k x , k y , k) via the following filtering function: Sensors 2024, 24, 2931 8 of 20 An important feature of the filtering function ( 12) is that it does not depend on lateral coordinates (x s , y s ) of scatterers and their scattering strengths K s , which is due to the factorized form of Equation (9).It also worth noting that for avoiding the occasional division by nearly zero values, in the denominator of Equation ( 12) some regularization term α > 0 should be added Rigorously speaking, the filtering function in Equation ( 12) depends on the wavenumber k = k n of every component of the illuminating optical field.Equation ( 12) also indicates that the filtering depends on the scatterer depth z = z s .However, we recall that unlike coordinate z s of s-th scatterer, which is a continuous quantity, the depth coordinate z q in the reconstructed OCT scans (see Equation ( 7)) is discreet.Thus, for transformation of OCT scans it makes sense to define the depth-dependent filtering function (12) for discreet depth points z q .Now, one can obtain the corrected (target) A-scan A (T) x 0 , y 0 , z q by transforming the initial A-scan A x 0 , y 0 , z q given by Equation ( 7) to the K-space, then multiplying the so-found spectrum by the filter function (12) and then again applying the Fourier transform to return to the corrected spectrum to the spatial domain.This procedure looks as follows: Now, one may recall that the spectral composition S(k n ) of the illuminating beam for OCT devices is well localized around the central wavenumber k 0 ≈ k n .Taking this into account, Equation (13) can be significantly simplified by eliminating the inverse and forward Fourier transforms over the wavenumbers k n .The latter can be substituted by the central wavenumber k 0 .Then, the simplified version of Equation (13) takes the form: In comparison with Equation (13), the simplified Equation ( 14) strongly reduces the computational expenses with only an insignificant influence on the accuracy of the entire correction procedure.
Thus, the derived Equations ( 13) and ( 14) (together with Equation ( 12) for the filtering function in the k-space) describe the OCT-scan transformation to the desired form.The transformation steps can be summarized in the diagram shown in Figure 3.In the next section, we will present some examples of application of the described transformation procedure.
Sensors 2024, 24, x FOR PEER REVIEW 9 of 22 form to return to the corrected spectrum to the spatial domain.This procedure looks as follows: ( ) Now, one may recall that the spectral composition of the illuminating beam for OCT devices is well localized around the central wavenumber . Taking this into account, Equation (13) can be significantly simplified by eliminating the inverse and forward Fourier transforms over the wavenumbers n k .The latter can be substituted by the central wavenumber 0 k .Then, the simplified version of Equation ( 13) takes the form: In comparison with Equation ( 13), the simplified Equation ( 14) strongly reduces the computational expenses with only an insignificant influence on the accuracy of the entire correction procedure.
Thus, the derived Equations ( 13) and ( 14) (together with Equation ( 12) for the filtering function in the k-space) describe the OCT-scan transformation to the desired form.The transformation steps can be summarized in the diagram shown in Figure 3.In the next section, we will present some examples of application of the described transformation procedure.

Conventional Refocusing
The first example demonstrating the efficiency of the above-derived relationships is their application for obtaining a deeper insight into the widely discussed problem of digital refocusing of OCT data.For digital shifting, the focus depth using a 3D set of OCT data obtained with a Gaussian beam, the following procedure based on several spectral transformations is known [25,26]: where ( ) , , 0 is the initial field and ( ) is the optical field after the numerical shift of the focal plane at a distance z Δ .In paper [25], this procedure was proposed using

Examples of Transformation
Procedures A(x 0 ,y 0 ,z q )→A (T) (x 0 ,y 0 ,z q ) 5.1.Conventional Refocusing The first example demonstrating the efficiency of the above-derived relationships is their application for obtaining a deeper insight into the widely discussed problem of digital refocusing of OCT data.For digital shifting, the focus depth using a 3D set of OCT data obtained with a Gaussian beam, the following procedure based on several spectral transformations is known [25,26]: where A 0 (x, y, z) is the initial field and A ∆z (x, y, z) is the optical field after the numerical shift of the focal plane at a distance ∆z.In paper [25], this procedure was proposed using heuristic arguments without rigorous derivation.In what follows, we clarify in detail the background approximations and demonstrate how this expression can be rigorously derived from the general Equation ( 12) for the filtering function and Equations ( 13) and ( 14) for spectral transformations based on this function.
Let us consider the illuminating beam U beam of a Gaussian form using the notations from work [11] assuming the focal depth z = z 0 : Here, W 0 is the radius of the focal waist at z = z 0 and the current beam width W(z, z 0 ) as a function of z is given by where λ = 2π/k.Assuming that the illuminating/receiving aperture is located at z = 0, W(0, z 0 ) ≡ W L is the initial radius of the beam at z = 0. Function Φ(z, ) describes the phase difference between the Gaussian beam and an ideally plane wave propagating the same distance (z − z 0 ).This phase term does not affect the beam focusing because it is independent of lateral coordinates, so in the following consideration it may be omitted.Quantity R(z, z 0 ) in Equation ( 16) is the current phase-front curvature so that R(0, z 0 ) is the initial curvature at the aperture.At this point, one may recall that the necessity of refocusing is important for highly focused OCT systems, for which W 0 ∼ λ near the focus and π • W 2 0 /(λz 0 ) << 1, so that Similarly, in view of π • W 2 0 /(λz 0 ) << 1, the relationship (17) between the focal radius W 0 and the beam radius W(0, z 0 ) ≡ W L at the aperture is also simplified Thus, omitting pre-factors (including the term with Φ(0, z 0 )) independent of lateral coordinates and, therefore, insignificant for focusing, the initial field distribution over the aperture at z = 0 can be written as Here, pre-factor C comprises independent lateral coordinate terms that are insignificant for focusing.The first exponential factor describes the distribution of the illuminating beam amplitude over the aperture plane z = 0, whereas the position of the focus corresponding to z = z 0 is determined by the wave-front curvature.For the Gaussian distribution (21) at the aperture, in the approximation that W L >> (λz 0 ) 1/2 (i.e., with the requirement that there are many Fresnel zones over the aperture), the angular spectrum g(k x , k y , k) has the form Now spectrum (22) can be used to find the filtering function (12).To this end, the auxiliary integral I(z s ; k, k x , k y ) given by Equation ( 10) should be first evaluated for the initial position z 0 of the illuminating-beam focus, and then target integral I (T) (z s , k, k x , k y ) should be similarly evaluated for any desired shifted focus distance z ′ 0 = z 0 − ∆z.In the general case, integral (10) can be evaluated numerically by taking the propagator function h s (. ..) in the full form given by Equation ( 2) in combination with function g(. ..) given by Equation (22).However, only for comparison with previously known results, we intentionally perform this evaluation analytically using the propagator h s (. ..) in the paraxial approximation instead of full Equation (2): In this approximation, the integrand in Equation ( 10) contains the following expression: For further consideration, we note that the last factor in Equation ( 24) independent of k x,y can be omitted.Both the angular spectrum g(k x , k y , k) given by Equation ( 22) and the propagator function h s (k x , k y , k) in Equation ( 24) contain phase factors with similar quadratic dependences on the lateral wave-vectors k x,y : exp i • z 0 (k , for which integral (10) can be found analytically.The resultant expression obtained from Equation (10) for integral I(. ..) (i.e., for the initial focus distance z = z 0 − z ′ o and scatterer depth z s ) contains the following factors In the target integral I (T) (. ..) we keep the same scatterer position z s (that is inherited from Equation (23) for h s (. ..)), but require that the new position z ′ 0 of the focal plane is shifted by ∆z = z 0 − z ′ 0 and assume that the beam radius and the aperture W ′ L is not necessarily equal to the initial W L .Thus, I (T) (. ..) takes the form functionally very similar to Equation ( 25) Finally, the filtering function (12), equal to the ratio I (T) (. ..)/I(. ..) takes the form It is clear that the refocusing factor in the earlier introduced Equation ( 15) differs from Equation (27) in which, besides the first phase exponential factor, an additional second exponential factor is present.This second exponential factor in Equation ( 27) disappears if Then, Equation ( 27) reduces to the same phase-only correction (in the spectral domain) as in Equation (15).However, this actually means that in the spatial domain, Equation (15) corresponds to the variable beam radius described by Equation (28).Due to this correction of the beam radius dependent on the refocusing depth, the focus waist after refocusing, based on Equation (15), keeps exactly the same radius as in the physical focus.This is clear from the relationship (20) between W 0 and W L .In contrast, if the beam radius in Equation ( 27) is invariable at the aperture, W L = W ′ L , then after refocusing, the focus-waist radius is either smaller (for z ′ 0 < z 0 ) or larger (for z ′ 0 > z 0 ) than in the physical focus at z = z 0 .We emphasize that the general Equation ( 12) for the filtering function is reduced to analytical Equation ( 27) and to previously known Equation (15), assuming the paraxial approximation and under the assumption that the curvature radius at the aperture is approximated as Note finally that, by analogy with general Equation ( 14), with good accuracy Equation ( 15) can also be represented in the simplified form This simplified form requires a much smaller number of calculations than Equation ( 15), and is especially computationally efficient.
Figure 4a demonstrates a simulated OCT B-scan representing a vertical chain of seven point-like scatterers imaged by a highly focused Gaussian beam, for which only the central scatterer located at the focus depth z s = z 0 = 252 µm is imaged with the maximal resolution.In this example, a central wavelength of 0.85 µm is chosen as quite a typical value for OCT devices (however, for another desirable wavelength, the results may be easily recalculated).Details of the usage of model Equations ( 5)- (7) for scan generation can be found in [19].Notice that these equations also readily allow one to generate OCT scans with a high density of scatterers.Examples of the latter can be found in [19], but in such scans, regions near and far from the focus depth visually do not clearly differ.Therefore, in Figure 4 and subsequent examples, we intentionally use only several distinct scatterers to clearly show the effects of beam transformations.Figure 4a clearly demonstrates that out of focus, the point-spread function becomes broader and gets a bit curved.The lateral size of the simulated scan in Figure 4 is chosen to be ~3-4 times greater than the illuminating-beam diameter in the out-of-focus regions to reduce the influence of boundaries on the results of the digital Fourier transformations required for refocusing.A further increase in simulated region size technically is not a problem, but this does not make much sense because in the too large field of view, the sizes of the initial and refocused scatterer images would be poorly seen for clear comparison.
Figure 4b shows the result of digital refocusing based on Equation ( 15) which actually performs depth-dependent adjusting of the optical-beam radius at the illuminating aperture and thus enables the depth-independent focus radius.Figure 4c,d show the results of refocusing with the use of the proposed filtering function (12) and Equation ( 14) for constructing refocused images under two assumptions: first, assuming the invariable beam width W L at the illuminating aperture (Figure 4c), and then invariable focus radius W 0 , but correspondingly corrected W L for the current focus depth according to Equation (20) (Figure 4d).
curved.The lateral size of the simulated scan in Figure 4 is chosen to be ~3-4 times greater than the illuminating-beam diameter in the out-of-focus regions to reduce the influence of boundaries on the results of the digital Fourier transformations required for refocusing.A further increase in simulated region size technically is not a problem, but this does not make much sense because in the too large field of view, the sizes of the initial and refocused scatterer images would be poorly seen for clear comparison.15) that looks as a phase-only correction but implicitly implies variable width of the beam at the aperture, such that . (c) is the result of refocusing based on the filtering function (12) and Equation ( 14) in which the focus depth is shifted, but the beam radius at the aperture is kept invariable , so that the focus radius  12) and ( 14), but requiring that  15) that looks as a phase-only correction but implicitly implies variable width of the beam at the aperture, such that is the result of refocusing based on the filtering function (12) and Equation ( 14) in which the focus depth is shifted, but the beam radius at the aperture is kept invariable W ′ L = W L0 , so that the focus radius is the result again based on Equations ( 12) and ( 14), but requiring that W ′ 0 = W 0 , such that condition W ′ L /z ′ 0 = W L /z 0 was used when calculating the filtering function (12).The central wavelength is λ = 0.85 µm, the initial focus depth is z 0 = 252 µm, and the radius of physical focus is W 0 = 1.9 µm.
To conclude this discussion, it can be mentioned that in Equations ( 15) and ( 29) the shift ∆z of the focus looks like a fixed value.However, to obtain refocusing of image A ∆z (x, y, z) at all depths z = z p (where p = 1 . . .N), one should recall that for pixels with different coordinate z, the focus shift ∆z from the initial position z 0 should be different, ∆z = z 0 − z.Bearing this condition in mind, for simultaneous refocusing of pixels at all depths z, Equation ( 29) can be rewritten as We point out that the matrix expression (30) performs a one-step transformation of the initial 3D set of pixelated OCT data A 0 (x, y, z) into another 3D array A ∆z (x, y, z).As a result, pixels at all depths z = z p become refocused without the need of additional operations (such as auxiliary transforms between the spatial and z-frequency domains, etc.).Actually, the refocused images in Figure 4 were obtained in such a way.

A Non-Gaussian and Non-Paraxial Example of Using the Spectral Filtering Function
In Section 5.1, we performed a comparison of OCT-image refocusing based on an earlier known correcting function of the phases of angular spectral components given by Equation ( 15), with the derived generalized spectral filtering function (12).For a clearer comparison, we evaluated the latter function using the conventional paraxial approximation and Gaussian beams, the properties of which are widely known.These approximations made it possible to analytically derive approximate forms of the auxiliary integrals (10) for the initial and target beams and thus reveal the connection between Equations ( 10)-( 12) and the earlier known results.
However, we emphasize that numerically performed transformations based on the filtering function given by Equations ( 10)-( 12) require neither paraxial approximation nor Gaussian beams.To illustrate the general nature of Equations ( 10)- (12), in this Section, we will use the full non-paraxial form of the propagator function (2) and intentionally demonstrate a bit of an exotic example of obtaining an OCT image for a non-diverging Bessel illuminating beam starting from an initial image acquired by a highly focused non-Gaussian beam.Namely, we start with a beam having a Lorentzian amplitude profile and spherical phase distribution over the aperture, for which amplitudes of spectral components k x,y decrease significantly slower than for a Gaussian beam, namely, ∝ exp(− k x,y ) instead of ∝ exp(−k 2 x,y ).For a Bessel beam, its amplitude profile at the aperture should be Gaussian, but the phase distribution should correspond to passing through an axicon lens [14,15] spectral filtering function (12) for such initial and target profiles is found digitally, using the numerical evaluation of auxiliary integrals (10) for the Lorentzian and Bessel-beam profiles.Figure 5 corresponds to the same set of scatterers as in Figure 4 and shows the vertical and horizontal OCT images corresponding to the initial Lorentzian beam and resultant Bessel beam.The intensities of scatterer images in Figure 5(a-1) are intentionally equalized to better show the out-of-focus scatterers.It is known that Bessel beams have noticeable sidelobes, but the central maximum well keeps its size independent of the depth and is comparable with the focal spot of the initial strongly focused beam.For clarity, this example is simulated in the absence of noise.The influence of noise on the results of such transformations and the necessity of the corresponding regularization procedures will be discussed in the next section.
Sensors 2024, 24, x FOR PEER REVIEW 14 of 22 beam.The intensities of scatterer images in Figure 5a-1 are intentionally equalized to better show the out-of-focus scatterers.It is known that Bessel beams have noticeable sidelobes, but the central maximum well keeps its size independent of the depth and is comparable with the focal spot of the initial strongly focused beam.For clarity, this example is simulated in the absence of noise.The influence of noise on the results of such transformations and the necessity of the corresponding regularization procedures will be discussed in the next section.z closer to the image bottom.The insets show the corresponding horizontal profiles.

Super-Refocusing
In Section 5.1, we demonstrated the use of the derived spectral filtering function ( 12) for digital refocusing of an initial 3D set of OCT data.The analysis based on Equation (12) revealed that, for ensuring depth-independent lateral resolution in the refocused OCT data, the same as in the physical focus, it was necessary to manipulate both phase curvature and the beam width at the aperture.This combined phase-amplitude manipulation is not obvious from the form of filtering function given by Equation (15), which looks like a phase-only correction in the spectral domain.The use of the proposed generalized filtering function (12)  Now, like in Section 5.1, we consider the most interesting for practice Gaussian beams and demonstrate how the filtering function (12) can be used for phase-amplitude manipulation with the initial OCT data to enable "super-refocusing".Possibilities of increasing resolution in optics in general have been analyzed for a long time (e.g., review [27]).Specifically for OCT, these issues were recently discussed in [28] using both theoretical argumentation based on the information capacity of optical images (similarly to [27]) and even some experimental examples of resolution enhancement.Such "super-refocusing" does not mean that the diffraction limit is overcome, nevertheless, proper

Super-Refocusing
In Section 5.1, we demonstrated the use of the derived spectral filtering function (12) for digital refocusing of an initial 3D set of OCT data.The analysis based on Equation (12) revealed that, for ensuring depth-independent lateral resolution in the refocused OCT data, the same as in the physical focus, it was necessary to manipulate both phase curvature and the beam width at the aperture.This combined phase-amplitude manipulation is not obvious from the form of filtering function given by Equation ( 15), which looks like a phase-only correction in the spectral domain.The use of the proposed generalized filtering function (12) in Section 5.1 made it possible to demonstrate that the choice of equal initial and target beam radii W L = W ′ L at the aperture resulted in a variable focal radius W ′ 0 , so that for z ′ 0 < z 0 , the W ′ 0 is somewhat smaller than W 0 in the physical focus.Now, like in Section 5.1, we consider the most interesting for practice Gaussian beams and demonstrate how the filtering function (12) can be used for phase-amplitude manipulation with the initial OCT data to enable "super-refocusing".Possibilities of increasing resolution in optics in general have been analyzed for a long time (e.g., review [27]).Specifically for OCT, these issues were recently discussed in [28] using both theoretical argumentation based on the information capacity of optical images (similarly to [27]) and even some experimental examples of resolution enhancement.Such "super-refocusing" does not mean that the diffraction limit is overcome, nevertheless, proper modification of the phase-amplitude beam profile can enable a significantly smaller target focus radius than in the physical focus, W T 0 < W 0 .In what follows, we use the derived generalized filtering function (12) to obtain such resolution enhancement and demonstrate the limiting role of measurement noises on the efficiency of such super-refocusing.
Recalling the relationship W 0 ≈ λz 0 /(πW L ) following from Equation ( 17) and the filtering function (27) derived for conventional refocusing of Gaussian beams in Section 5.1, it is clear that to enable invariable W T 0 < W 0 , one should choose in Equation ( 27) invariable ratio z ′ 0 /W ′ L < z 0 /W L for every refocusing depth z ′ 0 , and simultaneously the beam width W ′ L should be chosen broader than the initial beam radius.An example of the application of such procedures for super-refocusing is presented in Figure 6.In this example, the simulation is made in the presence of noise in the initial image, assuming a typical OCT SNR~30 dB (Figure 6(a-1,a-2)).The initial Gaussian illuminating beam is assumed to be weakly focused (with a focus radius of W 0 = 12.6 µm and the central wavelength λ = 0.85 µm is the same as in the previous examples).Now, instead of a vertical chain of individual scatterers, we consider a similar chain of pairs of point-like scatterers laterally separated by 12 µm.In the image formed by such a weakly focused beam, the scatterers in the pairs are not resolved (see the vertical B-scan in Figure 6(a-1) and the en-face image of one of the pairs at z = 480 µm in the inset).To illustrate this statement, Figure 7    The demonstrated increase in resolution is attained due to the effective increase in the illuminating beam width at the aperture, or, in other words, due to applying an amplitude mask, which causes flattening of the beam profile.However, flattening of the beam is performed at the expense of suppression of its amplitude near the beam axis.This effect certainly reduces the amplitude of the retained signals from the scatterers and decreases the resultant SNR.Nevertheless, as demonstrated in Figure 6(a-1,a-2), for quite a realistic initial (see Figure 6(c-1,c-2)).Thus, such super-refocusing can be applied to significantly improve the imaging resolution by re-processing already acquired datasets made with low-resolution OCT setups if the initial SNR is not too bad.
In terms of the angular spectrum of the beam, this improvement in resolution is due to digitally performed enhancing of the out-of-center angular spectral components.These side spectral components initially have smaller amplitudes and, therefore, in real measurement, are fairly noisy.Consequently, enhancement of these components may significantly increase the noise in the super-refocused image.Thus, to prevent unacceptable noise increases, some regularization procedures should be used when performing the signal transformation.Although quite sophisticated forms of regularization are used in some works (e.g., [29]), their detailed discussion is far beyond the scope of this paper.In view of this, we limit ourselves to the simple regularization method mentioned after the expression (12) for the spectral filter.Namely, to avoid division by near-zero values, we add a small regularization parameter α > 0 to the denominator in Equation ( 12).This parameter strongly affects the resultant noise level in the transformed OCT image.
To illustrate this statement, Figure 7 shows the simulated super-refocused profiles of the OCT signal through the pair of scatterers located at z = 480 µm, similarly to Figure 6(c-2), but using various values of α.The upper two rows Figure 7(a-1-a-5,b-1-b-5) show the simulated results in the absence of any measurement noises.All profiles are normalized to the peak amplitudes corresponding to the smallest α = 10 −6 , these amplitudes are taken as unity value (0 dB in Figure 7(a-1)).Figure 7(a-1-a-5)) demonstrate that a gradual increase in α by 6 orders of magnitude leads to a gradual decrease in the maxima (up to a 10 dB decrease for α = 1 in Figure 7(a-5), in comparison with the initial 0 dB for α = 10 −6 in Figure 7(a-1)).At the same time, the minimum in the center between the scatterers becomes less deep.Simultaneously, noise-like sidelobes around the scatterer images arise with increasing α.
Sensors 2024, 24, x FOR PEER REVIEW 17 of 22 scatterers becomes shallower.Figures 6 and 7 show that the improvement of resolution by a factor of 4 in this example leads to ~2-2.5 times a decrease in SNR for the optimal regulation parameter 2 10 − = α .Thus, the improvement of resolution after super-refocusing is obtained at the expense of decrease in SNR in the resultant image.If the initial SNR is not too small, after duly chosen regularization degree corresponding to the initial noise, the resultant SNR and the resolution of the image after super-refocusing may be still acceptable as in Figure 7(c-4,d-4).Although such super-refocusing improves the initial lateral resolution, it cannot overcome the diffraction limit, even if the initial image has a very high SNR.Furthermore, the digital enhancement of the effective numerical aperture of the illuminating beam may be limited by the physical lens aperture; the latter also limits the attainable effective value of the numerical aperture and, thus, the possible degree of super-refocusing.

Correction of Aberrations at the Illuminating/Receiving Aperture
Another problem for which the developed K-space framework can be applied is the correction of aberrations at the illuminating/receiving aperture.Here, we consider the The next two rows, Figure 7(c1-c5,d1-d5), show how the super-refocused images and their profiles look for the same set of α-values in the presence of moderate noise, the same as in the image before refocusing in Figure 6(a-1,a-2).It is clear that for too small α = 10 −6 and α = 10 −4 , the super-refocused images of scatterers in Figure 7(c-1,c-2,d-1,d-2) as in Figure 4a.It is clear that conventional refocusing based on Equation (15) (the middle row in Figure 8) makes the images of scatterers noticeably closer to the image in the physical focus.However, the so-refocused images still retain pronounced aberration-induced distortions of the point-spread function, as well as the lateral shift in the positions of scatterer images.The lowest row in Figure 8 (panels (a3)-(c3)) shows the results of refocusing combined with the compensation of the aberration function Q(x, y) = exp[−iAZ 1  3 (ρ, ϕ)] based on the spectral filtering described by Equations ( 12)-( 14).It is clear that the last row demonstrates both nearly perfect refocusing and compensation of the aberrations (the tiny residual distortions of <1% are related to the finite discretization step in the numerical evaluation of exact analytical expressions).15) (middle row) or refocusing combined with elimination of the aberrations using the spectral filter given by Equations ( 12) and ( 14

Discussion and Conclusions
In this paper, in the development of a series of earlier works related to the use of angular-spectrum approach to the description of OCT-image formation, we illustrate some possibilities of transformation of OCT data opened by this powerful approach.In the general sense, this approach develops the ideas of Fourier optics formulated several decades ago (see, e.g., textbooks [40,41]).For the realization of these principles, the progress in the development of powerful numerical packages (such as Matlab, Mathematica, etc.) opened rather convenient possibilities.
The main result of this paper is the derivation of a general form of the filtering function that enables various transformations of 3D sets of OCT data.The developed rather general and rigorous formulation of spectral-filtering procedures is summarized in Figure 3 and mathematically presented by Equations ( 10)- (14).To the best of our knowledge, they have not been derived before.The derived expressions in some special cases can be evaluated analytically.In more general cases, they are intended for numeri-  15) (middle row) or refocusing combined with elimination of the aberrations using the spectral filter given by Equations ( 12) and ( 14) (lower row).Column 1 shows the in-depth B-scans, column 2 shows the en face image of the scatter located at the depth of physical focus z = 252 µm and column 3 is the en face image of the scatterer located deeper at z = 480 µm.The small insets show the horizontal profiles of the scatterer images along the X-axis.

Discussion and Conclusions
In this paper, in the development of a series of earlier works related to the use of angular-spectrum approach to the description of OCT-image formation, we illustrate some possibilities of transformation of OCT data opened by this powerful approach.In the general sense, this approach develops the ideas of Fourier optics formulated several decades ago (see, e.g., textbooks [40,41]).For the realization of these principles, the progress in the development of powerful numerical packages (such as Matlab, Mathematica, etc.) opened rather convenient possibilities.
The main result of this paper is the derivation of a general form of the filtering function that enables various transformations of 3D sets of OCT data.The developed rather general and rigorous formulation of spectral-filtering procedures is summarized in Figure 3 and mathematically presented by Equations ( 10)- (14).To the best of our knowledge, they have not been derived before.The derived expressions in some special cases can be evaluated analytically.In more general cases, they are intended for numerical evaluation, for the realization of which such conventional assumptions as the paraxial approximation or the necessity of Gaussian beams are not required.The formulated approach, in particular, made it possible to clearly reveal conditions implicitly assumed in the conventionally discussed refocusing (compare the filtering functions in Equations ( 15) and ( 27)).We demonstrated that although this refocusing looks like the application of a phase-correcting factor in the spectral domain, actually this procedure corresponds to the digital correction of both phase and amplitude to shape the illuminating beam.Without such a simultaneous correction, the lateral size of the digitally refocused point-spread functions (i.e., the lateral resolution after refocusing) becomes dependent on the depth and may be either somewhat greater or smaller than in the initial focus (see Figure 4).To demonstrate the generality of the derived spectral filtering function ( 12), we presented an example of the transformation of a highly focused beam with an initial Lorentzian profile into a narrow Bessel beam.
Next, besides conventional refocusing, we demonstrated that the filtering function ( 12) allows one to simultaneously manipulate the phase and amplitude shape of the illuminating beam for "super-refocusing".Certainly, this does not yet mean overcoming the diffraction limit, but the lateral resolution after such super-refocusing may become several times higher than in the initial physical focus.The ultimate resolution attainable using such super-refocusing, similarly to the physical-focus radius, is limited by the well-known diffraction limit.In real-life situations, the attainable resolution is additionally affected by the level of noise, which may be the main limiting factor in super-refocusing rather than the diffraction limit.However, for quite a realistic OCT initial signal-to-noise ratio ~25-30 dB, it is still possible to obtain an increase in the lateral resolution ~several times, as illustrated in Figures 6 and 7.
Generally, the same procedures of digital filtering in the K-space may be used for other transformations of OCT signals, in particular, for compensation of aberrations as discussed in Section 5.3.Concerning other possibilities for analysis of OCT-signals and studies of various distorting factors, it was pointed out in [19,21] that the K-space description of OCT-scan formation makes it possible to simulate dispersion effects.
The procedures for forming OCT scans described here and in [19] enable a flexible and computationally efficient framework for realistic simulations of OCT signals in various situations, for which direct experimental tests are too difficult/expensive.All examples demonstrated in the manuscript required a few seconds of computation for the initial generation and transformations.Even for rather heavy simulations of forward problems, such as 3D data volumes containing large amounts of scatterers, computations can be made reasonably rapidly, although the required time increases proportionally with the number of scatterers.For example, with a volume of 128 × 128 × 128 pixels, corresponding to 192 × 192 microns laterally and 512 microns in depth with 250,000 scatterers (which imitate the density of biological cells ~5-6 microns in size), the simulations required 59 min using 12 cores of the AMD Ryzen 9 3900X CPU.However, for subsequent transformations of such an image, the computation time does not depend on the number of scatterers and requires only a few seconds, like in the above-presented examples with several scatterers.Numerous examples of problems in which such simulations can be used are mentioned in [19,21].In the latter work, specific examples related to studying non-Gaussian beams or such modalities as OCT-elastography and OCT-angiography are given.Overall, the K-space-based approach opens rather diverse possibilities for numerical modeling of OCT data and their digital transformations, as well as for deep analysis of OCT-signal features for usage in biomedical imaging and other applications.We also emphasize that the utilization of such simulated (but fairly realistic) OCT data opens unique possibilities for

Figure 1 .
Figure 1.Schematically shown geometry of the illuminating beam (not necessarily focused like in the figure) with the lateral coordinates ) , ( 0 0 y x of the axis directed to the tissue depth along the

Figure 1 .
Figure 1.Schematically shown geometry of the illuminating beam (not necessarily focused like in the figure) with the lateral coordinates (x 0 , y 0 ) of the axis directed to the tissue depth along the z-axis of the coordinate system.Coordinates (x s , y s , z s ) characterize positions of scatterers.
(1).Thus, the spectrum transformation during the back propagation is described by the same total spectral factor (3) for the forward propagation.Consequently, the transition from the spectral representation to the received signal amplitude is again given by the Fourier transform, similar to Equation (4).

Figure 2 .
Figure 2. Schematic elucidation of the fact that the forward propagation of the illuminating light and back-scattered signal propagation can be described in a symmetrical manner leading to Equation (6).Notice that such an equivalent scheme should be represented for every particular scatterer depth s z .

Figure 2 .
Figure 2. Schematic elucidation of the fact that the forward propagation of the illuminating light and back-scattered signal propagation can be described in a symmetrical manner leading to Equation (6).Notice that such an equivalent scheme should be represented for every particular scatterer depth z s .

Figure 3 .
Figure 3. Diagram showing main steps of transformation of the complex valued OCT data obtained for the initial field distribution ) ; , ( n L k y x U over the illuminating/receiving aperture to the desired form corresponding to a target field distribution

Figure 3 .
Figure 3. Diagram showing main steps of transformation of the complex valued OCT data obtained for the initial field distribution U L (x, y; k n ) over the illuminating/receiving aperture to the desired form corresponding to a target field distribution U (T) L (x, y; k n ) over the aperture.

Figure 4 .
Figure 4. Comparison of differently performed refocusing of a B-scan for a vertical chain of sub-resolution scatterers.Panel (a) is the initial B-scan with a highly focused beam simulated using Equation (7) for each single A-scan.(b) is the result of refocusing based on Equation (15) that looks as a phase-only correction but implicitly implies variable width of the beam at the aperture, such that is the result again based on Equations (

Figure 4 .
Figure 4. Comparison of differently performed refocusing of a B-scan for a vertical chain of subresolution scatterers.Panel (a) is the initial B-scan with a highly focused beam simulated using Equation (7) for each single A-scan.(b) is the result of refocusing based on Equation (15) that looks as a phase-only correction but implicitly implies variable width of the beam at the aperture, such thatW ′ L /z ′ 0 = W L /z 0 , which yields W ′ 0 = W 0 .(c)is the result of refocusing based on the filtering function(12) and Equation (14) in which the focus depth is shifted, but the beam radius at the aperture is kept invariable W ′ L = W L0 , so that the focus radiusW ′ 0 ̸ = const.(d)is the result again based on Equations (12) and (14), but requiring that W ′ 0 = W 0 , such that condition W ′ L /z ′ 0 = W L /z 0 was used when calculating the filtering function(12).The central wavelength is λ = 0.85 µm, the initial focus depth is z 0 = 252 µm, and the radius of physical focus is W 0 = 1.9 µm.

Figure 5 .
Figure 5. Example of transformation of non-Gaussian beams: initially strongly focused beam with a Lorentzian profile (row (a-1)-(a-3)) is transformed into a Bessel beam (row (b-1)-(b-3)).First column is for vertical B-scans, the other columns are horizontal scans through the depth 1 z of the

'
in Section 5.1 made it possible to demonstrate that the choice of equal initial and target beam radii W is somewhat smaller than 0 W in the physical focus.

Figure 5 .
Figure 5. Example of transformation of non-Gaussian beams: initially strongly focused beam with a Lorentzian profile (row (a-1-a-3)) is transformed into a Bessel beam (row (b-1-b-3)).First column is for vertical B-scans, the other columns are horizontal scans through the depth z 1 of the initial focus and z 2 closer to the image bottom.The insets show the corresponding horizontal profiles.

22 Figure 6 .
Figure 6.Comparison of conventional refocusing and super-refocusing for an initially weakly focused beam with 6 .12 0 = W μm for the same vertical positions of scatterers as in Figure 4.At each depth, a pair of scatterers laterally separated by 12 μm is located.Upper row shows color-coded B-scans with insets showing en-face images of scatterers located at z = 480 μm, for which the lateral profiles are shown in the lower raw.Panels (a-1,a-2) are not refocused and (b-1,b-2) show conventional refocusing similar to that in Figure 4. (c-1,c-2) show the result of super-refocusing with 4-fold increase in the lateral resolution, so that individual scatterers become clearly resolved.The noise level is especially clearly seen in the lower row.
shows the simulated super-refocused profiles of the OCT signal through the pair of scatterers located at 480 = z μm, similarly to Figure 6(c-2), but using various values of α .The upper two rows Figure 7(a-1-a-5,b-1-b-5)show the simulated results in the absence of any measurement noises.All profiles are normalized to the peak amplitudes corresponding to the smallest 6 10 − = α , these amplitudes are taken as unity value (0 dB in Figure7(a-1)).
Figure7(a-1-a-5)) demonstrate that a gradual increase in α by 6 orders of magnitude leads to a gradual decrease in the maxima (up to a 10 dB decrease for 1 a-1)).At the same time, the minimum in the center between the scatterers becomes less deep.Simultaneously, noise-like sidelobes around the scatterer images arise with increasing α .

Figure 6 .
Figure 6.Comparison of conventional refocusing and super-refocusing for an initially weakly focused beam with W 0 = 12.6 µm for the same vertical positions of scatterers as in Figure 4.At each depth, a pair of scatterers laterally separated by 12 µm is located.Upper row shows color-coded B-scans with insets showing en-face images of scatterers located at z = 480 µm, for which the lateral profiles are shown in the lower raw.Panels (a-1,a-2) are not refocused and (b-1,b-2) show conventional refocusing similar to that in Figure 4. (c-1,c-2) show the result of super-refocusing with 4-fold increase in the lateral resolution, so that individual scatterers become clearly resolved.The noise level is especially clearly seen in the lower row.

Figure 7 .
Figure 7. Illustration of how the regularization parameter influences the noise level in the super-refocused image with the same parameters as in Figure 5. Rows (a-1-a-5,b-1-b-5) show the "digital" noise caused by the regularization itself in the absence of other noises.Rows(c-1-c-5,d-1-d-5) correspond to the same moderate initial noise as in Figure5and demonstrate that the resultant noise after refocusing strongly depends on the regularization parameter.

Figure 7 .
Figure 7. Illustration of how the regularization parameter influences the noise level in the superrefocused image with the same parameters as in Figure 5. Rows (a-1-a-5,b-1-b-5) show the "digital" noise caused by the regularization itself in the absence of other noises.Rows (c-1-c-5,d-1-d-5) correspond to the same moderate initial noise as in Figure5and demonstrate that the resultant noise after refocusing strongly depends on the regularization parameter.

Sensors 2024 , 22 Figure 8 .
Figure 8. Illustration of distortions produced by the aberration in the form of Zernike polynomial ) , ( 1 3 ϕ ρ Z (upper row) and results of the image transformation using either conventional refocusing by Equation (15) (middle row) or refocusing combined with elimination of the aberrations using the spectral filter given by Equations (12) and (14) (lower row).Column 1 shows the in-depth B-scans, column 2 shows the en face image of the scatter located at the depth of physical focus 252 = z μm and column 3 is the en face image of the scatterer located deeper at ) (lower row).Column 1 shows the in-depth B-scans, column 2 shows the en face image of the scatter located at the depth of physical focus 252 = z μm and column 3 is the en face image of the scatterer located deeper at 480 = z μm.The small insets show the horizontal profiles of the scatterer images along the X-axis.

Figure 8 .
Figure 8. Illustration of distortions produced by the aberration in the form of Zernike polynomial Z 1 3 (ρ, ϕ) (upper row) and results of the image transformation using either conventional refocusing by Equation (15) (middle row) or refocusing combined with elimination of the aberrations using the spectral filter given by Equations (12) and (14) (lower row).Column 1 shows the in-depth B-scans, column 2 shows the en face image of the scatter located at the depth of physical focus z = 252 µm and column 3 is the en face image of the scatterer located deeper at z = 480 µm.The small insets show the horizontal profiles of the scatterer images along the X-axis.